Домашнее задание 3¶
0. Введение¶
Проблемой короновируса COVID-19 буквально сейчас заняты многие ученые по всему миру. В стороне не осталось и комьюнити Data scientist-ов, которые пытаются не просто предсказать развитие пандемии, а помочь определить наиболее значимые факторы, влияющие на распространение заразы. Например на небезызвестном Kaggle.
В данной домашней работе Вы будете использовать ежедневно обновляемые данные из репозитория, а конкретно time-series-19-covid-combined.csv.
Вам предстоит изучить имеющиеся данные, попробовать найти некоторые закономерности и ответить на вопросы. Ниже представлены задания как исследовательского характера, где приветствуется инициативность, так и чисто тренировочного харатера - где необходимо использовать изученные на семинаре инструменты для работы с временными рядами. Не забывайте сопровождать ваш код развернутыми комментариями и выводами. Чем их больше и чем они качественнее, тем лучше и выше вероятность получения максимальной оценки за домашнее задание. Не болейте!
In [1]:
!wget -O covid.csv https://raw.githubusercontent.com/datasets/covid-19/master/data/time-series-19-covid-combined.csv --no-check-certificate -q
In [1]:
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
pd.options.display.float_format = '{:.2f}'.format
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (15, 15)
In [200]:
df = pd.read_csv('covid.csv', parse_dates=['Date'])
df.head()
Out[200]:
1. Exploratory data analysis¶
1.1 Определите границы исследуемого периода, количество заболевших, погибших, выздоровивших. Отобразите на карте распространение болезни.¶
Границы исследования
In [201]:
df['Date'].agg(['min', 'max'])
Out[201]:
Количество заболевших, выздоровивших и погибших
In [202]:
df[df['Date'] == df['Date'].max()][['Confirmed', 'Recovered', 'Deaths']].sum()
Out[202]:
Карта
Данные логарифмированы, иначе будет одна яркая точка -- США. На карте доступна анимация, а так же масштабирвоание.
In [5]:
import folium
import folium.plugins as plugins
In [6]:
data = []
__df = df[(~df['Confirmed'].isna()) & (df['Province/State'] != 'Diamond Princess')].copy()
__df['Confirmed'] = np.log(1 + __df['Confirmed'])
__df['Confirmed'] = __df['Confirmed'] / __df['Confirmed'].max()
__df = __df[['Date', 'Lat', 'Long', 'Confirmed']].groupby('Date').agg(lambda x: list(x))
for date, row in __df.iterrows():
row = list(zip(*row))
row = [[x[0], x[1], np.clip(x[2], 0.00001, 1)] for x in row if x[2] > 0]
data.append(row)
m = folium.Map([0., 0.], zoom_start=2)
hm = plugins.HeatMapWithTime(data,
index=list(__df.index.strftime('%Y-%m-%d').values),
auto_play=True,
min_speed=1,
)
hm.add_to(m)
m
Out[6]:
1.2 Определите количество стран, которые подтвердили случаи заражения. Визуализируйте рост количества стран, подтвердивших случаи заражения по дням (или неделям).¶
Примечание. В анализе временных рядов иногда необходимо изменить частоту дискретизации (Resampling). Такая операция может понадобится, когда имеющихся данных недостаточно или их, наоборот, слишком много. Может не устраивать имеющаяся частота или просто хочется посмотреть на данные с более общего ракурса.
Различают два вида изменения частоты: повышение (Upsampling) и понижение (Downsampling). При повышении временной ряд пересчитывается с низкой частоы на более высокую частоту (например от годовой до месячной частоты). В таком случае процесс включает в себя заполнение или интерполяцию появившихся пропусков в данных. При понижении временной ряд передискретизируется с высокой частоты на низкую (наример с еженедельной на месячную частоту). Это включает в себя агрегацию существующих данных.
In [203]:
data = df[df['Confirmed'] > 0].groupby('Date')['Country/Region'].nunique()
data.name = 'Approved countries'
data.index = data.index.strftime('%Y-%m-%d')
data.plot.bar(figsize=(15, 7), title='Количество зараженных стран по дням');
In [8]:
data = df[df['Confirmed'] > 0].groupby(pd.Grouper(key='Date', freq='W-MON'))['Country/Region'].nunique()
data.name = 'Approved countries'
data.index = data.index.strftime('%Y-%m-%d')
data.plot.bar(figsize=(15, 7), title='Количество зараженных стран по неделям');
1.3 Визуализируйте мировые тренды по числу заболевших | погибших | выздоровивших. Постройте графики, которые отражают летальность от вируса в разных странах или на разных континентах.¶
In [9]:
data = df.groupby('Date')[['Confirmed', 'Deaths', 'Recovered']].sum()
data.index = data.index.strftime('%Y-%m-%d')
data.plot(figsize=(15, 7), title='Количество заболевших | погибших | выздоровивших');
In [10]:
__df = df[(df['Date'] == df['Date'].max()) & (df['Confirmed'] > 0)].copy()
__df = __df.groupby('Country/Region')[['Deaths', 'Confirmed']].sum()
__df['cfr'] = __df['Deaths'] / __df['Confirmed']
__df.sort_values('cfr', inplace=True)
__df = __df[__df['Confirmed'] > 100]
# __df = __df[__df['cfr'] > __df['cfr'].mean()]
__df['cfr'].plot.barh(figsize=(7, 30), title='Летальность в разных странах');
1.4 Визуализируйте и сравните, как распространялась болезнь в разных странах или континентах по дням с первого зарегистрированного случая заболевания. Определите количество дней для каждой страны со дня первого зарегестрированного случая.¶
Подчистим данные
In [205]:
df_c = df.groupby(['Date', 'Country/Region'])[['Confirmed', 'Recovered', 'Deaths']].sum().reset_index()
df_c
Out[205]:
In [206]:
region_start = df_c[df_c['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start.name = 'd_start'
__df = pd.merge(df_c[df_c['Confirmed'] > 0], region_start, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start']).dt.days
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df
Out[206]:
In [207]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if i < 10:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число зараженных')
plt.title('Распространение болезни в разных странах с первого подтвержденного случая');
In [208]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Active', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Active']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Active_max', ascending=False)
).iterrows()):
if i < 10:
color = plt.plot(row['n_day'], row['Active'], label=f'{i+1:>2}.{name:<15}- total {row["Active_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['Active'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Active'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число активных больных')
plt.title('Число активно больных в разных странах с первого подтвержденного случая');
In [209]:
(
__df
.groupby('Country/Region')['n_day']
.max()
.sort_values(ascending=False)
.plot
.barh(figsize=(7, 30), title='Сколько дней прошло с первого заражения')
);
Построим (для интереса) аналогичные графики от 100 зараженного¶
In [16]:
region_start_100 = df_c[df_c['Confirmed'] >= 100].groupby('Country/Region')['Date'].min()
region_start_100.name = 'd_start100'
__df = pd.merge(df_c[df_c['Confirmed'] >= 100], region_start_100, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start100']).dt.days
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df
Out[16]:
In [18]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if i < 10:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel('Число зараженных')
plt.title('Распространение болезни в разных странах с 100 подтвержденного случая');
In [19]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Active', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Active']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Active_max', ascending=False)
).iterrows()):
if i < 10:
color = plt.plot(row['n_day'], row['Active'], label=f'{i+1:>2}.{name:<15}- total {row["Active_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['Active'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Active'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число активных больных')
plt.title('Число активно больных в разных странах с первого подтвержденного случая');
In [20]:
(
__df
.groupby('Country/Region')['n_day']
.max()
.sort_values(ascending=False)
.plot
.barh(figsize=(7, 30), title='Сколько дней прошло с 100 заражения')
);
Особой разницы между графиками не заметно, но от 100 зараженного выглядит более гладким
Если смотреть на Китай и Германию, то на переход от начала заражения до пика проходит меньше времени, чем от пика до спада.
1.5 Визуализируйте и оцените как развивалась ситуация в России по сравнению с остальными миром.¶
In [21]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if 1 <= i < 15:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel('Число зараженных')
plt.ylim([-4000, 1.1 * __df.groupby('Country/Region')['Confirmed'].max().sort_values()[-2]])
plt.title('Распространение болезни в разных странах с 100 подтвержденного случая');
plt.tight_layout()
Данные на 2020-04-23¶
Количество подтвержденных случаев в России близко к аналогичным в Бразилии, Канаде, Китаю и Ирану.
Если сравнивать темпы роста, то Россия похожа на UK, Францию, Италию, Германию.
По предыдущим графикам, у России еще все впереди.
1.6 Определите топ-10 стран, которые лидируют по имеющимся в датасете показателям. Сравните насколько быстро развивается ситуация в этих странах, определите корреляцию.¶
In [22]:
for col in ['Confirmed', 'Recovered', 'Deaths', 'Active']:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day'], row[col], label=f'{i+1:>2}.{name:<15}- total {row[f"{col}_max"]:>8,d}' if i < 20 else None)[0].get_color()
if 0 <= i < 10:
plt.scatter(row['n_day'][-1], row[col][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row[col][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel(f'Число {col}')
# plt.ylim([-4000, 1.1 * __df.groupby('Country/Region')[col].max().sort_values()[-2]])
plt.title(f'Распространение болезни ({col}) в разных странах с 100 подтвержденного случая');
plt.tight_layout()
plt.show()
Данные на 2020-04-23¶
Все показатели, кроме выздоровивших похожи. Лидирует США. Во некоторых странах ситуация намного хуже чем в Китае.
По количеству выздоровивших неожиданно лидирует Германия.
In [23]:
import seaborn as sns
def plot_corr(D, size):
corr = D.corr()
# corr = np.abs(corr)
f, ax = plt.subplots(figsize=(size, size))
cmap = plt.cm.Oranges
sns.heatmap(corr, cmap=cmap,
xticklabels=corr.columns,
yticklabels=corr.columns)
from IPython.display import display
In [24]:
for col in ['Confirmed', 'Recovered', 'Deaths', 'Active']:
print(f'Корреляция по признаку {col}')
display(
(__df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
.head(10)
)
.T
.apply(lambda x: pd.Series(x.loc[col]))
.corr()
.style
.background_gradient(cmap='coolwarm')
.set_precision(2)
)
Данные на 2020-04-23¶
По всем показателям, кроме Recovered Китай не похож на остальные страны. В топе стран высокая кореляция по наблюдаемым величинам. По Recovered все страны похожи
1.7 Для первых 5 стран из топ-10 и России постойте сравнительные графики. Изучите как изменялись значения отношений погибших/выздоровивших, погибших/заболевших, заболевших/общее число жителей, ваш вариант...
In [25]:
!wget -O population.csv https://raw.githubusercontent.com/datasets/population/master/data/population.csv --no-check-certificate -q
In [26]:
top_c = list(set(df_c.groupby('Country/Region')['Confirmed'].max().sort_values(ascending=False).head(5).index) | {'Russia'})
top_c
Out[26]:
In [27]:
population = pd.read_csv('population.csv', sep=',')
population['Country Name'].replace({'United States': 'US', 'Russian Federation': 'Russia', 'Iran, Islamic Rep.': 'Iran'}, inplace=True)
population = population[population['Country Name'].isin(top_c)]
population = population.loc[population.groupby('Country Name')['Year'].idxmax().values].set_index('Country Name')['Value']
population.name = 'population'
population.index.name = 'Country/Region'
set(population.index) - set(top_c)
Out[27]:
In [28]:
__df = df_c[df_c['Country/Region'].isin(top_c)]
# region_start_100 = df_c[df_c['Confirmed'] >= 100].groupby('Country/Region')['Date'].min()
# region_start_100.name = 'd_start100'
__df = __df.join(region_start_100, on='Country/Region').join(population, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start100']).dt.days
__df = __df[__df['n_day'] > 0]
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df['Deaths/Recovered'] = __df['Deaths'] / (__df['Recovered'] + 1) # prevent divide by 0
__df['Deaths/Confirmed'] = __df['Deaths'] / __df['Confirmed']
__df['Confirmed/population'] = __df['Confirmed'] / __df['population']
__df.head()
Out[28]:
In [29]:
for col in ['Deaths/Recovered', 'Deaths/Confirmed', 'Confirmed/population', 'Active']:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row[col], label=f'{i+1:>2}.{name:<10}')[0].get_color()
plt.scatter(row['n_day'][-1], row[col][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row[col][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel(f'Число {col}')
plt.title(f'Распространение болезни ({col}) в разных странах с 100 подтвержденного случая');
plt.tight_layout()
plt.show()
Данные на 2020-04-23¶
Deaths/Recovered: С начала заболевания много людей погибает, а количество выздоровивших растет с некоторой задержкой, поэтому есть всплеск данного показателя в начале болезни в странах.
Deaths/Confirmed: В России подозрительно низкий показатель смертности. Если построить аналогичный график для 25 стран, то у России данный показатель также низок, на уровне Саудовской Аравии (в отличии от России, у нее он падал со временем). По всей видимости, в России по другому считают стасистику погибших (умерших от вируса, а не умерших с вирусом). Также, возможно, занижают показатели для более низкого значения смертности
Confirmed/population: Видно, что в США хоть и больше всего зараженных, но население страны так же большое. По данному показателю лидирует Испания. Также заметно, что заболело в каждой стране менее 1% населения, а в России менее 0,1%.
Active: В Германии заболевание идет на спад, в Италии, Испании и Франции рост прекращается. По России трудно говорить о начале пика заболевших. В США все плохо.
1.8 Выделите временной ряд по одному из целевых признаков. Для выделенного временного ряда:
- Представьте временной ряд в виде последовательности разницы между значениями. Рассчитайте лаги, выведите их на графике.
- Попробуйте декомпозировать временной ряд. Постройте график тренда.
- Что вы можете сказать о стационарности временного ряда?
In [30]:
__df = pd.concat([df_c[['Date', 'Country/Region', 'Confirmed']], df_c.shift(df_c.groupby('Date').size().unique()[0])['Confirmed']], axis=1).dropna()
__df.columns = ['Date', 'Country/Region', 'Confirmed', 'Confirmed_lag']
__df['lag'] = __df['Confirmed'] - __df['Confirmed_lag']
region_start = __df[__df['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start.name = 'd_start'
__df = pd.merge(__df[__df['Confirmed'] > 0], region_start, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start']).dt.days
__df['day_of_week'] = __df['Date'].dt.dayofweek
__df
Out[30]:
In [31]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['lag', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['lag']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('lag_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day'], row['lag'], label=f'{i+1:>2}.{name:<15}- total {row["lag_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['lag'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['lag'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Прирост больных')
plt.title('Число прироста новых случаев в разных странах');
Ряд нестационарен, заметны некоторые циклические колебания (по дням недели?).
Теоретически не может быть стационарным, когда-то все выздоровят. В Китае количество новых случаев давно стало нулевым
1.9 Для любого интересующего Вас временного ряда постройте графики ACF и PACF. Сделайте выводы.
In [32]:
from statsmodels.tsa.stattools import acf, pacf
In [33]:
pd.Series(acf(__df[__df['Country/Region'] == 'Russia']['lag'])).plot.bar(figsize=(14,7), title='ACF');
In [34]:
pd.Series(pacf(__df[__df['Country/Region'] == 'Russia']['lag'])).plot.bar(figsize=(14,7), title='PACF');
оптимальный параметры для модели ARMA (1, 2)
1.10. Изучаемый Вами датасет отражает текущую мировую статистику распространения болезни. Однако эти данные не дают нам возможноси предположить о возможных факторах, которые влюяют на скорость распространения, летальности, заразность и тд. Поэтому в наш датасет можно добавить новые признаки из других наборов данных, которые могут не только помочь при прогнозировании временного ряда, но и определить наиболее значимые факторы.¶
- Добавьте новый признак "количество дней со дня первого случая заражения". Добавьте еще не менее 5 новых признаков в датасет, которые могут быть значимым фактором распространения вируса. Сохраните полученный датасет, вы можете его использовать во второй части домашнего задания.
- Исследуйте как добавленные признаки соотносятся с трендами распространения короновируса в некоторых странах. Сделайте выводы.
Примеры признаков - количество дней со начала карантина, конинент, плотность населения, количество употребляемого алкоголя, летучих мышей, алкоголя, крокодилов на душу населения, количество туристов в год, уровень безработицы, температура, количество игроков в Plague Inc., средний возраст или продолжительность жизни и так далее. Вы можете использовать любой найденный датасет. Главное, чтобы у Вас получилось скомбинировать данные.
Примеры датасеты:
In [35]:
population = pd.read_csv('population.csv', sep=',')
population['Country Name'].replace(
{
'United States': 'US',
'Russian Federation': 'Russia',
'Iran, Islamic Rep.': 'Iran',
'Czech Republic': 'Czechia',
'Egypt, Arab Rep.': 'Egypt',
'Korea, Rep.': 'Korea, South',
'Slovak Republic': 'Slovakia',
}, inplace=True)
population = population.loc[population.groupby('Country Name')['Year'].idxmax().values].set_index('Country Name')['Value']
population.name = 'population'
population.index.name = 'Country/Region'
In [210]:
__df = df_c.copy()
__df['cfr'] = __df['Deaths'] / __df['Confirmed']
__df = pd.concat([__df, __df.shift(df_c.groupby('Date').size().unique()[0])['cfr']], axis=1).dropna()
cols = list(__df.columns)
cols[-1] = cols[-1] + '_lag'
__df.columns = cols
__df['lag'] = __df['cfr'] - __df['cfr_lag']
region_start1 = __df[__df['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start1.name = 'd_start1'
region_start100 = __df[__df['Confirmed'] > 100].groupby('Country/Region')['Date'].min()
region_start100.name = 'd_start100'
region_start1000 = __df[__df['Confirmed'] > 1000].groupby('Country/Region')['Date'].min()
region_start1000.name = 'd_start1000'
__df = (
__df[__df['Confirmed'] > 0]
.join(region_start1, on='Country/Region', how='inner')
.join(region_start100, on='Country/Region', how='inner')
.join(region_start1000, on='Country/Region', how='inner')
.join(population, on='Country/Region', how='inner')
)
__df['active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df['n_day1'] = (__df['Date'] - __df['d_start1']).dt.days
__df['n_day100'] = (__df['Date'] - __df['d_start100']).dt.days
__df['n_day1000'] = (__df['Date'] - __df['d_start1000']).dt.days
__df['confirmed_rate'] = __df['Confirmed'] / __df['population']
__df['recovered_rate'] = __df['Recovered'] / __df['population']
__df['deaths_rate'] = __df['Deaths'] / __df['population']
__df['active_rate'] = __df['active'] / __df['population']
__df['day_of_week'] = __df['Date'].dt.dayofweek
__df = __df[__df['n_day100'] >= 0]
__df.drop(columns=['cfr_lag', 'd_start1', 'd_start100', 'd_start1000', 'population'], inplace=True)
__df.reset_index(drop=True, inplace=True)
__df.to_csv('covid_features.csv')
__df.to_pickle('covid_features.pkl')
__df
Out[210]:
In [149]:
plt.figure(figsize=(20, 20), facecolor='w')
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
plot_num = 1
for i, (country, active) in enumerate(__df.groupby('Country/Region')['active'].max().sort_values(ascending=False).head(10).iteritems()):
___df = __df[__df['Country/Region'] == country]
for param in ['n_day1', 'n_day100', 'n_day1000', 'confirmed_rate', 'recovered_rate', 'deaths_rate', 'active_rate', 'day_of_week']:
plt.subplot(10, 8, plot_num)
if i == 0:
plt.title(param, size=18)
if param == 'n_day1':
plt.ylabel(country, fontsize=18)
if ___df.dtypes[param] == np.int64:
___df.groupby(param)['lag'].mean().plot.bar()
else:
___df.groupby(param)['lag'].mean().plot()
plt.xticks(())
plt.yticks(())
plot_num += 1
На графиках представлена зависимость прироста смертности от параметра дня в разных странах
1.11 Выводы¶
выводы были после каждого пункта
2. Прогнозирование временного ряда¶
Теперь вам предстоит построить несколько моделей и спрогнозировать временной ряд для такого показателя как летальности от вируса. Летальность можно рассчитывать по формуле:
$$ CFR = \frac{Deaths}{Confirmed} $$Для каждого типа модели сделайте несколько прогнозов для не менее 3 разных стран/регионов/другое (на ваш обоснованный выбор). Главное, чтобы каждый из типов моделей прознозировал одинаковый набор данных, чтобы в конце можно было сделать выводы о качестве работы той или иной модели.
Метрика качества RMSE. Не забудьте разбить данные временной ряд на данные для обучения и тестовые.
Не забывайте сопровождать ваш код комментариями, графиками и выводами.
In [150]:
import datetime
from sklearn.metrics import mean_squared_error
In [151]:
def rmse(y_true, y_pred):
return np.sqrt(mean_squared_error(y_true, y_pred))
In [152]:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df[__df['n_day1000'] > 0]
.groupby('Country/Region')[['lag', 'n_day1000']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['n_day1000']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('n_day1000_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day1000'], row['lag'], label=f'{i+1:>2}.{name:<15}')[0].get_color()
plt.scatter(row['n_day1000'][-1], row['lag'][-1], color=color, alpha=1)
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 1000 зараженного')
plt.ylabel('Прирост летальности')
plt.title('Число прироста летальности в разных странах');
In [153]:
split_date = __df['Date'].max() - datetime.timedelta(10)
print(split_date)
train = __df[__df['Date'] < split_date].copy()
test = __df[__df['Date'] >= split_date].copy()
In [42]:
all_countries = list(train['Country/Region'].unique())
top_countries = ['US', 'Russia', 'Spain']
In [43]:
def measure_quality(df, name):
per_country = (
df
.groupby('Country/Region')
.agg(lambda x: list(x))
.apply(lambda x: rmse(x['cfr'], x['cfr_r']) if np.isnan(x['cfr_r']).sum() == 0 else np.nan, axis=1)
)
per_country.name = 'rmse'
world = pd.merge(per_country, population, left_index=True, right_index=True, how='inner')
per_country['World'] = (world['rmse'] * world['population']).sum() / world['population'].sum()
per_country.name = name
return per_country
In [44]:
class Experimens:
def __init__(self, top_cols):
self._df = pd.DataFrame(columns=['World'] + top_cols)
def add(self, data, name):
quality = measure_quality(data, name)
self._df = self._df.append(quality, verify_integrity=True)
return self.show
@property
def show(self):
return self._df.copy()
In [45]:
ex = Experimens(top_countries)
2.1 Прогнозирование средним (Mean Constant Model)¶
In [46]:
def plot_predict(predict, title=None):
plt.figure(figsize=(14, 7))
for country in top_countries:
_tr = train[train['Country/Region'] == country]
_te = predict[predict['Country/Region'] == country]
color = plt.plot(_tr['Date'], _tr['cfr'], label=country + ' train')[0].get_color()
plt.plot(_te['Date'], _te['cfr'], label=country + ' test', color=color, marker='+')
plt.plot(_te['Date'], _te['cfr_r'], label=country + ' predict', color=color, marker='o')
plt.title(title)
plt.legend();
In [47]:
predict = train.groupby('Country/Region')['cfr'].mean()
predict = test[['Date', 'Country/Region', 'lag', 'cfr']].join(predict, on='Country/Region', rsuffix='_r')
plot_predict(predict, 'Константное предсказание')
ex.add(predict, 'const')
Out[47]:
Mean const lag¶
In [48]:
def flatten(df):
res = []
for index, row in df.iterrows():
for r2 in zip(*row):
new_row = {
df.index.name: index
}
for keys, values in zip(row.keys(), r2):
new_row[keys] = values
res.append(new_row)
return pd.DataFrame(res)
def predict_to_normal(data, data_start):
data_start = data_start[data_start['Date'] == data_start['Date'].max()].set_index('Country/Region')[['cfr']]
df = data[['Date', 'Country/Region', 'lag', 'lag_r']].groupby('Country/Region').agg(lambda x: list(x))
df['cfr'] = df['lag'].map(lambda x: list(np.array(x).cumsum()))
df['cfr_r'] = df['lag_r'].map(lambda x: list(np.array(x).cumsum()))
df.drop(columns=['lag', 'lag_r'], inplace=True)
df = flatten(df).join(data_start, rsuffix='_begin', on='Country/Region')
df['cfr'] = df['cfr'] + df['cfr_begin']
df['cfr_r'] = df['cfr_r'] + df['cfr_begin']
df.drop(columns=['cfr_begin'], inplace=True)
return df.copy()
In [49]:
predict = train.groupby('Country/Region')['lag'].mean()
predict = test[['Date', 'Country/Region', 'lag']].join(predict, on='Country/Region', rsuffix='_r')
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Константный лаг')
ex.add(predict, 'const lag')
Out[49]:
2.2 Линейная модель (Linear Trend Model)¶
In [50]:
import statsmodels.formula.api as smf
from tqdm import tqdm_notebook as tqdm
In [51]:
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = smf.ols('lag ~ n_day1', data=train[train['Country/Region'] == country]).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.predict(predict[index])
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Линейный лаг')
ex.add(predict, 'linear')
Out[51]:
2.3 Экспоненциальное сглаживание (Exponential smoothing)¶
Попробуйте так же применить Double exponential smoothing или Triple exponential smoothing.
In [52]:
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing
In [53]:
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = SimpleExpSmoothing(train[train['Country/Region'] == country].reset_index()['lag']).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.forecast(sum(index)).values
predict = predict_to_normal(predict, train)
plot_predict(predict, 'SimpleExpSmoothing')
ex.add(predict, 'SimpleExpSmoothing')
Out[53]:
In [54]:
from IPython.display import clear_output
In [55]:
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = ExponentialSmoothing(train[train['Country/Region'] == country].reset_index()['lag'], trend='additive').fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.forecast(sum(index)).values
clear_output()
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Double ExpSmoothing')
ex.add(predict, 'Double ExpSmoothing')
Out[55]:
2.4 Скользящее среднее (MA - Moving average)¶
Поэксперементируйте с гиперпараметрами модели.
In [56]:
from statsmodels.tsa.arima_model import ARIMA
In [57]:
for param in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (0, 0, param))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA(0, 0, {param}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA(0, 0, {param})')
ex.add(predict, f'ARIMA(0, 0, {param})');
ex.show
Out[57]:
2.5 Авторегрессионная модель (AR - Autoregressive model)¶
Поэксперементируйте с гиперпараметрами модели.
In [58]:
for param in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (param, 0, 0))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({param}, 0, 0) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({param}, 0, 0)')
ex.add(predict, f'ARIMA({param}, 0, 0)');
ex.show
Out[58]:
2.6 ARMA¶
Поэксперементируйте с гиперпараметрами модели.
In [59]:
for _AR in range(1, 4):
for _MA in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (_AR, 0, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({_AR}, 0, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({_AR}, 0, {_MA})')
ex.add(predict, f'ARIMA({_AR}, 0, {_MA})');
ex.show
Out[59]:
2.7 Интегрированная модель авторегрессии — скользящего среднего (Auto Regressive Integrated Moving Average)¶
In [60]:
_AR = 1
_MA = 1
for _I in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (_AR, _I, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({_AR}, {_I}, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({_AR}, {_I}, {_MA})')
ex.add(predict, f'ARIMA({_AR}, {_I}, {_MA})');
ex.show
Out[60]:
2.8 SARIMA¶
In [61]:
from statsmodels.tsa.statespace.sarimax import SARIMAX
In [62]:
_AR = 1
_I = 1
_MA = 1
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = SARIMAX(_tr, order=(_AR, _I, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index)).values
except:
bad += 1
if bad:
print(f'SARIMAX({_AR}, {_I}, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'SARIMAX({_AR}, {_I}, {_MA})')
ex.add(predict, f'SARIMAX({_AR}, {_I}, {_MA})')
Out[62]:
2.9 Библиотека Prophet¶
Попробуйте использовать библиотеку Prophet для предсказания временного ряда. Документация.
In [64]:
from fbprophet import Prophet
In [81]:
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
for country in tqdm(all_countries):
_tr = train[train['Country/Region'] == country].reset_index()[['Date', 'lag']]
_tr.columns = ['ds', 'y']
model = Prophet(weekly_seasonality=True)
result = model.fit(_tr)
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model.predict(model.make_future_dataframe(periods=sum(index)))['yhat'][-sum(index):].values
clear_output()
predict = predict_to_normal(predict, train)
plot_predict(predict, f'Phrofet')
ex.add(predict, f'Phrofet')
Out[81]:
2.10 Линейная модель с регрессией (Linear Regression Model)¶
Какие признаки оказались наиболее значимыми?
In [132]:
predict = test.copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = smf.ols('lag ~ n_day1 + n_day100 + n_day1000 + day_of_week', data=train[train['Country/Region'] == country]).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.predict(predict[index])
if country == 'Russia':
__res = model_linear.summary()
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Линейная модель с регрессией')
ex.add(predict, 'Линейная модель с регрессией')
Out[132]:
In [133]:
__res
Out[133]:
Наиболее значимые -- число дней после 100 и 1000 заражения, то есть когда ситуация по динамике стабилизировалась.
2.11 (Бонус) Используйте любую другую известную вам модель для предсказания.
In [135]:
from catboost import CatBoostRegressor
In [185]:
def df2normal(df, delay=7, col='lag', for_date=None):
df = df.sort_values('Date').groupby('Country/Region').agg(lambda x: list(x))
delay = datetime.timedelta(delay)
result = []
for country, __df in tqdm(df.iterrows(), total=len(df)):
for l_date, l_nday1, l_nday100, l_nday1000, l_lag in zip(__df['Date'], __df['n_day1'], __df['n_day100'], __df['n_day1000'], __df[col]):
if for_date and for_date != l_date:
continue
row = {
'date': (l_date - train['Date'].min()).days,
'dayofweek': l_date.dayofweek,
'n_day1': l_nday1,
'n_day100': l_nday100,
'n_day1000': l_nday1000,
'lag': l_lag,
'country': country,
}
for r_date, r_nday1, r_nday100, r_nday1000, r_lag in zip(__df['Date'], __df['n_day1'], __df['n_day100'], __df['n_day1000'], __df[col]):
if l_date - delay <= r_date < l_date:
delta = (l_date - r_date).days
row[f'date_{delta}'] = (r_date - train['Date'].min()).days
row[f'dayofweek_{delta}'] = r_date.dayofweek
row[f'n_day1_{delta}'] = r_nday1
row[f'n_day100_{delta}'] = r_nday100
row[f'n_day1000_{delta}'] = r_nday1000
row[f'lag_{delta}'] = r_lag
result.append(row)
result = pd.DataFrame(result)
return result
In [227]:
df3 = df2normal(__df, col='lag')
x_train = df3[df3['date'] < (split_date - train['Date'].min()).days].copy()
x_test = df3[df3['date'] >= (split_date - train['Date'].min()).days].copy()
y_train = x_train['lag']
y_test = x_test['lag']
x_train.drop(columns=['lag'], inplace=True)
x_test.drop(columns=['lag'], inplace=True)
In [228]:
clf_params = {
'depth': 6,
'num_trees': 10000,
'random_strength': 3,
'bagging_temperature': 1.02,
'cat_features': ['country'],
'eval_metric': 'RMSE',
'random_seed': 42,
'use_best_model': True,
'od_type': 'Iter',
'od_wait': 200,
'task_type': 'CPU',
}
In [229]:
clf = CatBoostRegressor(**clf_params)
clf.fit(x_train, y_train, eval_set=(x_test, y_test), plot=True, verbose=False)
Out[229]:
In [255]:
clf.get_feature_importance(prettified=True)
Out[255]:
In [230]:
df4 = __df.copy()
df4['lag'] = np.nan
for date in x_test['date'].unique():
date = train['Date'].min() + datetime.timedelta(int(date))
df4_ = df2normal(df4, col='lag', for_date=date)
df4_.drop(columns='lag', inplace=True)
df4.loc[df4['Date'] == date, 'lag'] = clf.predict(df4_)
In [232]:
predict = test.copy()
predict['lag_r'] = df4[~df4['lag'].isna()]['lag']
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Catboost')
ex.add(predict, 'Catboost')
Out[232]:
2.11 Выводы¶
In [234]:
ex.show.to_pickle('exp_results.pkl')
In [242]:
models = {k: 0 for k in ex.show.index}
top_models = {k: 0 for k in ex.show.index}
In [243]:
for col in ex.show.columns:
for i, model in enumerate(ex.show.sort_values(col).index):
if i == 0:
top_models[model] += 1
models[model] += i
Топ моделей по количеству первых мест:
In [250]:
pd.Series(top_models).sort_values(ascending=False)
Out[250]:
Рейтинг моделей по итоговой позиции:
In [252]:
pd.Series(models).sort_values(ascending=True)
Out[252]:
Как можно заметить, catboost победил все модели по качеству предсказания (выборка делилась по времени, данные для предсказания каждый раз обновлялись на основе предсказаний).
Тем не менее, некоторые алгоритмы близко подходят по качеству. Возможно, catboost выиграл за счет обобщающей способности по всем странам, возможно ряд не совсем классический для задач временных рядов.
Кроме того можно провести дополнительные эксперименты по предсказанию для catboost на cfr, а не лаги; попробовать применить LSTM и другие RNN.
Тем не менее, для каждой задачи возможно лучше будет использовать свой алгоритм, который лучше всего описывает данный ряд.